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which  depends  on  the  n-vector  state  x(t),  the  m-vector  control 
u(t),  and  the  p-vector  parameter  it  .\  The  state  is  given  at 
the  initial  point.  At  the  final  point,  the  state  and  the 
parameter  are  required  to  satisfy  q scalar  relations.  Along 
the  Interval  of  integration,  the  state,  the  control,  and  the 
parameter  are  required  to  satisfy  n scalar  differential 
equations.  Problem  P2  differs  from  Problem  PI  in  that  the 
state,  the  control,  and  the  parameter  are  required  to  satisfy 
k additional  scalar  relations  along  the  interval  of  integration. 
Algorithms  of  the  sequential  gradient-restoration  type 
are  given  for  both  Problem  PI  and  Problem  P2. 


Problem  P2  enlarges  dramatically  the  number  and  variety 
of  problems  of  optimal  control  which  can  be  treated  by 
gradient-restoration  algorithms,  i Indeed,  by  suitable  trans- 
formations, almost  every  known  problem  of  optimal  control  can 
be  brought  into  the  scheme  of  Problem  P2.  This  statement  applies,  for 
instance  to  the  following  situations:  (i).  problems 
with  control  equality  constraints,  (ii)  problems  with 
state  equality  constraints,  (iii)  problems  with  state- 
derivative  equality  constraints,  (iv)  problems  with  control 
inequal i ty  constra ints,  (v)  problems  with  state  inequality 
constraints,  and  (vi)  problems  with  state-derivative 
inequality  constraints. 

Eight  numerical  examples  are  presented  to  illustrate  the 
performance  of  the  algorithms  associated  with  Problem  PI 
and  Problem  P2.  The  numerical  results  show  the  feasibility 
as  well  as  the  convergence  characteristics  of  these 
algorithms. 
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1 . Introduction 

In  every  branch  of  science,  engineering,  and  economics, 
there  exist  systems  which  are  controllable,  that  is,  they  can 
be  made  to  behave  in  different  ways  depending  on  the  will  of 
the  operator.  Every  time  the  operator  of  a system  exerts  an 
option,  a choice  in  the  distribution  of  the  quantities  control- 
ling the  system,  he  produces  a change  in  the  distribution  of 
the  states  occupied  by  the  system  and,  hence,  a change  in  the 
final  state.  Therefore,  it  is  natural  to  pose  the  following 
question:  Among  all  the  admissible  options,  what  is  the  par- 
ticular option  which  renders  the  system  optimum?  As  an  example, 
what  is  the  option  which  minimizes  the  difference  between  the 
final  value  and  the  initial  value  of  an  arbitrarily  specified 
function  of  the  state  of  the  system?  The  body  of  knowledge 
covering  problems  of  this  type  is  called  calculus  of  variations 
or  optimal  control  theory.  As  stated  before,  applications 
occur  in  every  field  of  science,  engineering,  and  economics. 

It  must  be  noted  that  only  a minority  of  current  problems 
can  be  solved  by  purely  analytical  methods.  Hence,  it  is 
important  to  develop  numerical  techniques  enabling  one  to  solve 
optimal  control  problems  on  a digital  computer.  These  numeri- 
cal techniques  can  be  classified  into  two  groups:  first-order 
methods  and  second-order  methods.  First-order  methods  (or 
gradient  methods)  are  those  techniques  which  employ  at  most 
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the  first  derivatives  of  the  functions  under  consideration. 
Second-order  methods  (or  quasilinearization  methods)  are  those 
techniques  which  employ  at  most  the  second  derivatives  of  the 
functions  under  consideration. 

Both  gradient  methods  and  quasilinearization  methods  re- 
quire the  solution  of  a linear,  two-point  or  multi-point 
boundary-value  problem  at  every  iteration.  This  being  the 
case,  progress  in  the  area  of  numerical  methods  for  differen- 
tial equations  is  essential  to  the  efficient  solution  of 
optimal  control  problems  on  a digital  computer. 

In  this  paper,  we  review  recent  advances  in  the  area  of 
gradient  methods  for  optimal  control  problems.  Because  of 
space  limitations,  we  make  no  attempt  to  cover  every  possible 
technique  and  every  possible  approach,  a material  impossibility 
in  view  of  the  large  number  of  publications  available.  Thus, 
except  for  noting  the  early  work  performed  by  Kelley  (Refs.  1-2) 
and  Bryson  (Refs.  3-6) , we  devote  the  body  of  the  paper  to  a 
review  of  the  work  performed  in  recent  years  by  the  Aero- 
Astronautics  Group  of  Rice  University  (Refs.  7-34). 

Also  because  of  space  limitations,  we  treat  only  single- 
subarc problems.  More  specifically,  we  consider  two  classes 
of  optimal  control  problems,  called  Problem  PI  and  Problem  P2 
for  easy  identification. 

Problem  Pi  consists  of  minimizing  a functional  I which 
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depends  on  the  n-vector  state  x(t) , the  m-vector  control 
u(t),  and  the  p-vector  parameter  tt.  The  state  is  given  at 
the  initial  point.  At  the  final  point,  the  state  and 
the  parameter  are  required  to  satisfy  q scalar  relations. 

Along  the  interval  of  integration,  the  state,  the  control, 
and  the  parameter  are  required  to  satisfy  n scalar  differ- 
ential equations.  Problem  P2  differs  from  Problem  Pi  in 
that  the  state, the  control,  and  the  parameter  are  required  to 
satisfy  k additional  scalar  relations  along  the  interval  of 
integration.  Algorithms  of  the  sequential  gradient-restoration 
type  are  given  for  both  Problem  Pi  and  Problem  P2. 

1.1.  Approach.  The  approach  taken  is  a sequence  of  two- 
phase  cycles,  composed  of  a gradient  phase  and  a restoration 
phase.  The  gradient  phase  involves  one  iteration  and  is  de- 
signed to  decrease  the  value  of  the  functional,  while  the  con- 
straints are  satisfied  to  first  order.  The  restoration  phase 
involves  one  or  more  iterations,  and  is  designed  to  force 
constraint  satisfaction  to  a predetermined  accuracy,  while  the 
norm  squared  of  the  variations  of  the  control  and  the  parameter 
is  minimized,  subject  to  the  linearized  constraints. 

The  principal  property  of  the  algorithms  presented  here 
is  that  a sequence  of  feasible  suboptimal  solutions  is  produ- 
ced. In  other  words,  at  the  end  of  each  gradient-restoration 
cycle,  the  constraints  are  satisfied  to  a predetermined  accu- 
racy. Therefore,  the  values  of  the  functional  I corresponding 
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to  any  two  elements  of  the  sequence  are  comparable. 

The  stepsize  of  the  gradient  phase  is  determined  by  a 
one-dimensional  search  on  the  augmented  functional  J,  while 
the  stepsize  of  the  restoration  phase  is  obtained  by  a one- 
dimensional search  on  the  constraint  error  P.  The  gradient 
stepsize  and  the  restoration  stepsize  are  chosen  so  that  the 
restoration  phase  preserves  the  descent  property  of  the  gra- 
dient phase.  As  a consequence,  the  value  of  the  functional  I 
at  the  end  of  any  complete  gradient-restoration  cycle  is 
smaller  than  the  value  of  the  same  functional  at  the  beginning 
of  that  cycle. 

1.2.  Time  Normalization.  A time  normalization  is  used 
in  order  to  simplify  the  numerical  computations.  Specifically, 
the  actual  time  9 is  replaced  by  the  normalized  time  t = 9/x , 
which  is  defined  in  such  a way  that  t = 0 at  the  initial  point 
and  t = 1 at  the  final  point.  The  actual  final  time  t,  if  it  is 
free,  is  regarded  as  a component  of  the  vector  parameter  ir  to 
be  optimized.  In  this  way,  an  optimal  control  problem  with 
variable  final  time  is  converted  into  an  optimal  control  prob- 
lem with  fixed  final  time. 


5 


AP-16 


1.3.  Notation.  In  this  paper,  vector-matrix  notation  is 
used  for  conciseness. 

Let  t denote  the  independent  variable,  and  let  x(t),  u(t) , 
it  denote  the  dependent  variables.  The  time  t is  a scalar,  the 
state  x(t)  is  an  n-vector,  the  control  u(t)  is  an  m-vector, 
and  the  parameter  tt  is  a p-vector.  All  vectors  are  column 
vectors . 

Let  h(x,u,TT,t)  denote  a scalar  function  of  the  arguments 
x,u,tt  ,t.  The  symbol  h^  denotes  the  n-vector  function  whose 
components  are  the  partial  derivatives  of  the  function  h with 
respect  to  the  components  of  the  vector  x.  Analogous  defini- 
tions hold  for  h and  h . 

u tt 

Let  oj(x,u,7r,t)  denote  an  r-vector  function  of  the  arguments 
x,u,iT,t.The  symbol  denotes  the  nxr  matrix  function  whose 
elements  are  the  partial  derivatives  of  the  components  of 
the  vector  m with  respect  to  the  components  of  the  vector  x. 
Analogous  definitions  hold  for  the  symbols  w and 

The  dot  sign  denotes  derivative  with  respect  to  the  time, 
that  is,  x = dx/dt.  The  symbol  T denotes  transposition  of 
vector  or  matrix.  The  subscript  0 denotes  the  initial  point, 
and  the  subscript  1 denotes  the  final  point. 
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1.4.  Outline . Section  2 contains  the  statements  of 
Problem  Pi  and  Problem  P2.  Section  3 gives  a description  of 
the  sequential  gradient-restoration  algorithm.  Section  4 
discusses  the  determinations  of  the  basic  functions  for  the 
gradient  phase  and  the  restoration  phase.  Section  5 considers 
the  determination  of  the  stepsizes  for  the  gradient  phase 
and  the  restoration  phase.  A summary  of  the  sequential 
gradient-restoration  algorithm  is  presented  in  Section  6. 

The  experimental  conditions  are  given  in  Section  7.  The 
numerical  examples  for  Problem  PI  are  given  in  Section  8 ; and 
the  numerical  examples  for  Problem  P2  are  given  in  Section  9. 
Finally,  the  discussion  and  the  conclusions  are  presented  in 


Section  10. 
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2.  Statement  of  the  Problems 


Problem  Pi.  This  problem  consists  of  minimizing  the  functional 

r1 

1=  l f (x,u,ir,  t)  dt  + [g  (x,tt,  t)]  (1) 

J 0 

with  respect  to  the  state  x(t),  the  control  u(t),  and  the 
parameter  tt  which  satisfy  the  differential  constraints 

X - 0 (X,U,TT,t)  = 0,  0 < t < 1 , (2) 

the  initial  conditions 

x (0) = given,  (3) 

and  the  final  conditions 

[>J>  (x,ir,t)]  1 = 0.  (4) 

In  Eqs.  (l)-(4),  the  quantities  I,f,g  are  scalar,  the  function 
<t>  is  an  n-vector,  and  the  function  4>  is  a q-vector.  Eqs. 

(2) -(4)  constitute  the  feasibility  equations  for  Problem  Pi. 


Problem  P2.  This  problem  is  an  extension  of  Problem  PI, 


which  arises  because  of  the  inclusion  of  the  nondifferential 
constraints 


S (x , u , tt  , t ) = 0,  0 < t < 1,  (5) 

to  be  satisfied  everywhere  along  the  interval  of  integration. 
Here,  the  function  S is  a k-vector,  k < m.  Eqs.  (2) -(5)  cons- 
titute the  feasibility  equations  of  Problem  P2. 

Problem  P2  enlarges  dramatically  the  number  and  variety 
of  problems  of  optimal  control  which  can  be  treated  by 
gradient-restoration  algorithms.  Indeed,  by  suitable  trans- 
formations , almost  every  known  problem  of  optimal  control  can 
be  brought  into  the  scheme  of  Problem  P2.  This  statement 
applies,  for  instance,  to  the  following  situations:  (i)  prob- 
lems with  control  equality  constraints,  (ii)  problems  with 
state  equality  constraints, (iii) problems  with  state-derivative 
equality  constraints,  (iv)  problems  with  control  inequality 
constraints,  (v)  problems  with  state  inequality  constraints, 
and  (vi)  problems  with  state-derivative  inequality  constraints. 
For  an  illustration  of  the  scope  and  range  of  applicability 
of  Problem  P2,  the  reader  is  referred  to  Ref.  19  and  Refs. 

25-29. 

2.1.  Remark . For  both  Problem  Pi  and  Problem  P2,  the 
number  of  final  conditions  q must  satisfy  the  following  relation 
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q < n+  p*  < n + p , (6) 

where  the  symbol  p*  denotes  the  number  of  components  of 
the  parameter  tt  present  in  the  final  conditions. 

2.2.  Remark.  Problem  Pi  car.  be  regarded  as  a particular 
case  of  Problem  P2,  which  arises  by  deleting  Eq.  (5).  This 
being  the  case,  the  analytical  derivations  presented  here 
refer  only  to  Problem  P2.  The  corresponding  analytical  deri- 
vations for  Problem  Pi  can  be  obtained  by  setting 


o 

III 

cn 

(7) 

5 0, 

S = 0 

IT 

(8) 

in  the  equations  of  Problem  P2.  However,  the  differentiation 
between  Problems  PI  and  P2  is  invoked  later  on  in  the  paper, 
in  the  section  dealing  with  the  solution  of  the  linear,  two- 
point  boundary-value  problem  (LTP-BVP) . This  is  necessary 
for  computational  efficiency. 

2.3.  Augmented  Functional.  From  calculus  of  variations, 
it  can  be  seen  that  Problem  P2  is  one  of  the  Bolza  type,  which 

4 

can  be  recast  as  that  of  minimizing  the  augmented  functional 


I 
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J=  l [f  + Xi(x-'t>)+pis]dt+(g  + u •i1)^ 
0 


= (X 


•1 

? 1 T 

x)  Q + (f-X^t 
0 


pTS  - XTx)  dt  + (g  + pT'J;)  ^ , (9) 


subject  to  (2) -(5).  In  Eq.  (9),  X (t)  is  a variable  Lagrange 
multiplier  (an  n-ve.tor),  p(t)  is  a variable  Lagrange  multi- 
plier (a  k-vector) , and  p is  a constant  Lagrange  multiplier 
(a  q-vector) . 

2.4.  First-Order  Conditions.  Let  the  multipliers  X(t), 
p (t) , p be  chosen  consistently  with 


X-f  +<|>X-Sp  = 0, 

X X X 


0 < t < 1 , 


(X  + gx  + ^xp)  1 = 0. 

Then,  the  optimal  control  u(t)  and  parameter  tt  satisfy 

the  following  relations: 


f — 4>  X + S p = 0,  0 < t < 1 , 

u u u - - 


(10) 


(ii) 


(12] 


In  Eq.  (9) , it  is  tacitly  assumed  that  the  initial  conditions 
(3) are  satisfied.  The  second  form  of  Eq.  (9)  arises  after 
the  customary  integration  by  parts  is  performed. 


, 
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\ (f^  - ^X  + S^pjdt  + (g^  + iPttu)1  = 0.  (13) 

J 0 

< 

Eqs. (10) - (13) constitute  the  optimality  conditions  for  Problem  P2. 


t 
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2.5.  Two-Point  Boundary-Value  Problem.  The  system 

(2)-  (5)  and  (10) -(13)constitutes  a nonlinear,  two-point  boundary- 
value  problem  in  which  the  unknowns  are  the  functions  x(t), 
u(t),  it  and  the  multipliers  X (t)  , p(t),  u.  Only  for  particular 
cases,  closed-form  solutions  are  possible.  In  general,  nume- 
rical methods  must  be  employed. 

Depending  on  whether  these  numerical  methods  employ  at 
most  the  first  derivatives  or  at  most  the  second  derivatives 
of  the  functions  under  consideration,  two  classes  of  algorithms 
can  be  developed:  first-order  algorithms  (also  called 
gradient  methods)  and  second-order  algorithms  (also  called 
quasilinearization  methods) . As  stated  in  the  introduction, 
only  first-order  algorithms  are  considered  here. 

2.6.  Performance  Indexes.  When  solving  Problem  P2  on  a 
digital  computer,  it  is  necessary  to  define  convergence  in  a 
numerical  sense.  In  this  connection,  let  tt»e  norm  squared  of 
a vector  y be  defined  as 

N (y)  = yTy.  (14) 


/ 
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Let  P and  Q denote  the  scalar  performance  indexes 


P=\  N(x-<f)dt+|  N(S)dt  + N (^ ) x , 


(15) 


1 

0 


Q=|  N(X  - fx  + <)>XA  - Sxp)  dt  + ^ N (fu  - ^uX  + Sup)  dt 


+N 


(f7T_  V + S7rp)dt  + + 


+ N(X+gx+ ^xy)1,  (16! 


which  measure  the  errors  in  the  constraints  and  the  optimality 
conditions,  respectively.  Observe  that 
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3 . Sequential  Gradient-Restoration  Algorithm 

The  technique  employed  is  characterized  by  a sequence  of 
two-phase  cycles,  composed  of  a gradient  phase  and  a restor- 
ation phase.  The  gradient  phase  is  started  only  when  Ineq. 
(19-1)  is  satisfied;  it  involves  one  iteration  and  is  designed 
to  decrease  the  value  of  the  functional  I or  the  augmented 
functional  J,  while  the  constraints  are  satisfied  to  first 
order.  The  restoration  phase  is  started  only  when  Ineq. 

(19-1)  is  violated;  it  involves  one  or  more  iterations,  each 
designed  to  decrease  the  constraint  error  P,  while  the  norm 
squared  of  the  variations  of  the  control  u(t)  and  the  parame- 
ters n is  minimized.  The  restoration  phase  is  terminated 
whenever  the  constraints  are  satisfied  to  a predetermined  ac- 
curacy, that  is,  whenever  Ineq.  (19-1)  is  satisfied. 

A complete  gradient-restoration  cycle  is  designed  so  that 
the  value  of  the  functional  I decreases  while  the  constraints 
are  satisfied  to  the  accuracy  (19-1)  both  at  the  beginning  and 
at  the  end  of  the  cycle.  Finally,  the  algorithm  as  a whole  is 
terminated  whenever  Ineqs.  (19)  are  satisfied  simultaneously. 

3.1.  Notation.  For  any  iteration  of  the  gradient  phase 
or  the  restoration  phase,  the  following  terminology  is  adopted; 
x(t),  u(t),  7r  denote  the  nominal  functions;  x(t),  u(t),  n 
denote  the  varied  functions;  and  Ax(t),  Au(t),  Air  denote  the 
displacements  leading  from  the  nominal  functions  to  the  varied 
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functions.  These  quantities  satisfy  the  relations 

x(t)=x(t)+Ax(t),  u(t)=u(t)+Au(t),  tt  = tt  + Att.  (20) 

Let  a be  a positive  number  representing  the  stepsize 
(either  the  gradient  stepsize  or  the  restoration  stepsize) . 

Then,  we  define  the  displacements  per  unit  of  stepsize  as 
follows : 

A ( t)  = Ax  ( t)  /a  , B(t)  = Au(t)/a,  C=AiT/a.  (21) 

Upon  combining  (20)  and  (21),  we  see  that 

x(t)  = x(t)  + aA(t)  , u (t)  = u (t)  + aB  (t)  , tf  = tt  + aC.  (22) 


3.2.  Desired  Properties.  The  functions  Ax(t),  Au(t), 

Att  must  be  determined  so  as  to  produce  some  desirable  effect 
at  every  iteration,  namely,  the  decrease  of  the  functionals  I, 
and/or  J,  and/or  P.  Thus,  the  following  descent  properties 
are  required: 

I < I,  and/or  J < J,  and/or  P < P,  (23) 
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where  I,  J,  P are  associated  with  the  nominal  functions  and  I, 
J,  P are  associated  with  the  varied  functions.  In  turn,  the 
functions  A(t) , B(t),  C are  chosen  so  that 


51 < 0,  and/or  5J  < 0,  and/or  5P  < 0, 


where  the  symbol  5(...)  denotes  the  first  variation.  Then,  by 

choosing  the  stepsize  a sufficiently  small,  the  satisfaction 

of  relations  (23)  is  guaranteed.  Ineqs.  (23-1) , (23-2)  and 

(24-1),  (24-2)  characterize  the  gradient  phase,  while  Ineqs. 

(23-3)  and  (24-3)  characterize  the  restoration  phase. 

3.3.  First  Variations.  Next,  we  give  the  expressions 

for  the  first  variations  of  the  functionals  I,  J,  P;  after 

simple  manipulations,  omitted  for  the  sake  of  brevity,  they 
, 6,7 


take  the  form 


5 I/a 


■E 


rp  fp  rp  rn  rp 

(fxA+fuB+f,C)dt+  (*xA  + ^C)l' 


^Implicit  in  Eqs.  (25)-(27)  is  the  assumption  A(0)  =0. 

^The  first  variation  of  the  augmented  functional  J is  computed 
by  varying  the  functions  x(t),  u(t),  if,  while  holding  the 
multipliers  X(t),p(t),  y unchanged. 
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6J/a  = \ (-a  + f - <p  X + S p)  1 Adt  + 

XX  X 
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I (f  -<t>  X + s 

I U U 

J 0 


p) TBdt 


rr1 
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and 

+ 

1 <f„-V +V,dt+<'?„  + Vh 

.Jo 

ct L(A  + W»  \ 

, (26) 


6P/2a=\  (x  - 0)T(A  - (f^A  - - <J>^C)dt 


0 
*1 

T T T T I rp  rn 

+ \ S (S  A+S  B + S C)dt+  (^A+ 

1 X U 7T  X ’ 


(27: 


For  the  purposes  of  this  paper,  Eqs . (25) -(27)  must  be  completed 
by  the  following  relation: 


1 

0 


2 I T T 

K/ct  = \ B Bdt  + C C , 


(23: 


which  constitutes  a measure  of  the  overall  change  of  the  con- 
trol and  the  parameter. 

3.4.  Remark.  Clearly,  every  iteration  of  either  the 
gradient  phase  or  the  restoration  phase  includes  two  distinct 
operations:  (a)  the  determination  of  functions  A(t),  B(t),  C 
consistent  with  the  first  variation  requirements  (24);  and 
(b)  the  determination  of  the  stepsize  a consistent  with  the 
total  variation  requirements  (23)  . 


I 
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4 . Determination  of  the  Basic  Functions 

There  exist  an  infinite  number  of  combinations  of  fun- 
ctions A ( t),  B ( t),  C capable  of  satisfying  the  first-variation 
inequalities  (24),  subject  to  the  linearized  constraints.  In 
order  to  arrive  at  a unique  combination  of  functions,  some 
additional  requirement  must  be  imposed.  This  is  done  through 
the  formulation  of  the  following  auxiliary  minimization  prob- 
lems . 

Problem  p3.  For  the  gradient  phase,  minimize  the  linear 
functional  (25),  with  respect  to  the  perturbations  A(t),  B(t), 

C which  satisfy  the  linearized  constraints 


A - ^x  ~ = 0 , 0 < t < 1 , 

(29) 

sta  + stb  + sTc  = 0 , 0 < t < 1 , 

X U 7T  - 

(30) 

A ( 0 ) = 0, 

(31) 

oj£a  + ^c)x  = o, 

(32) 

and  the  quadratic  isoperimetric  constraint  (28). 

Note  the  absence  of  forcing  terms  from  Eqs.  (29)  — (32)  . 

This  implies  that  the  nominal  functions  characterizing  the 
gradient  phase  satisfy  the  constraints  (2)- (5)  within  the  pre- 
selected accuracy  (19-1) . 


Problem  P4.  For  the  restoration  phase,  minimize  the 
quadratic  functional  (28)  , with  respect  to  the  perturbations 
A ( t ),  B ( t),  C which  satisfy  the  linearized  constraints 


I 
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A - *TA  - 0TB  - '{>TC  + (x  - <|>)  = 0 , 

r X U IT 


0 < t < 1 , 


T T T 

S A + S B + S C + S = 0, 
x u rr 


0 < t < 1, 


A ( 0 ) = 0, 


(^A  + + *)1  = 0.  (3( 

Note  that  forcing  terms  are  absent  from  Eq.  (35)/  but  are 
present  in  Eqs.(33),  (34),  (36).  This  implies  that  the  nominal 
functions  characterizing  the  restoration  phase  satisfy  the 


initial  conditions  (3),  but 


violate  one  or  more  of  the 


remaining  constraints  (2),  (4),  (5) 


Indeed,  it  is 


the  purpose  of  the  restoration  phase  to  correct  these  viola- 
tions, while  causing  the  least  possible  disturbance  in  the 
system.  This  is  the  significance  of  the  least-square  criterion 
(28)  . 

4.1.  First-Order  Conditions.  Problems  P3  and  P4  are 
variational  problems  of  the  Bolza  type,  each  governed  by  a 
different  set  of  optimality  conditions. 

For  Problem  P3,  let  the  multipliers  A(t),  p(t),  p be  chosen 
consistently  with 


X " fx  + (pxX  " V = °' 


0 < t < 1, 


(\  + gx  + = 0. 
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Then,  the  control  perturbation  B(t)  and  parameter  perturbation 
C satisfy  the  following  relations: 


B+fu-V  + Sup  = 0' 


0 < t < 1 , 


c+  \ (fn  - <PnX  + s^pjdt  + (gu  + ^u)  ± = 0 


For  Problem  P4  , let  the  multipliers  X(t),  p(t),  u be  chosen 
consistently  with 


X + <t>X-Sp  = 0,  0 < t < 1 , 

xx  - 


(X  + '^xu ) ^ = o . 


Then,  the  control  perturbation  B(t)  and  parameter  perturbation 
C satisfy  the  following  relations: 


B - <?> u X + SuP  = 0, 


0 < t < 1, 


C+  I (-^A+S^pjdt  + (i^u)  x = 0.  (' 

J 0 

4.2.  Linear,  Two-Point  Boundary-Value  Problem.  For  the 
gradient  phase,  Problem  P3  is  governed  by  the  feasibility 
equations  (29)  — (32)  and  the  optimality  conditions  (37)  — (40) . 
For  the  restoration  phase,  Problem  P4  is  governed  by  the 
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feasibility  equations  (33) -(36)  and  the  optimality  conditions 
(41) -(44).  The  form  of  these  feasibility  equations  and 
optimality  conditions  is  such  that  the  systems  governing 
Problems  P3  and  P4  can  be  embedded  in  a single  linear  system. 
For  compactness,  as  well  as  to  facilitate  programming,  this 
point  of  view  is  taken  here,  and  the  single  system  governing 
both  the  gradient  phase  and  the  restoration  phase  is  written 
as  follows: 


A - 0TA  - 0TB  - 4>TC  + k (x  - <f> ) = 0 , 
x u TT  r 


STA  + S^B  + S^C  + k S = 0, 
X U 7T  r 


0 < t < 1 , 


0 < t < 1 , 


A ( 0 ) = 0, 


(<p*A+<p*C  + k’pL  = 0, 

X II  X.  -L 


X-kf  + 'p  X - S p = 0,  0<t<l, 

g x x x - 


(UVx  + V)l3°' 


B + kgfu‘ V + Sup  = 0'  °-  t-1' 


C+  (*g^"*7TX  + V)dt  + (*ggir  + Vi>1“°* 

J 0 

The  constants  k and  kr  appearing  in  (45)- (52)  take  the 
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following  values: 

gradient  phase, 

V1' 

k =0; 
r 

(53) 

restoration  phase. 

o 

II 

tp 

kr  = l. 

(54) 

For  given  nominal  functions  x(t),  u(t),  it  and  given  constants 
kg  and  k^,  Eqs.  (45) — (52)  define  the  functions  A(t),  B(t),  C 
and  the  multipliers  X ( t ),  p(t),  y.  As  can  be  seen,  we  are  in  the 
presence  of  a linear,  two-point  boundary-value  problem  (LTP-BVP), 
which  can  be  solved  independently  of  the  value  assigned  to  the 
stepsize  a. 

^ In  principle, 

the  LTP-BVP  (45) -(52)  can  be  discussed  simultaneously  for  both 
Problem  Pi  and  Problem  P2.  However,  for  computational  effi- 
ciency, it  is  better  to  separate  the  discussion  of  Problem  Pi 
from  that  of  Problem  P2.  This  is  because  the  LTP-BVP  for 
Problem  Pi  can  be  solved  executing  q + 1 independent  sweeps  of 
the  differential  system,  while  the  LTP-BVP  for  Problem  P2 
requires  the  execution  of  n + p + 1 independent  sweeps.  Here, 
q is  the  number  of  final  conditions,  n is  the  dimension  of  the 
state  vector,  and  p is  the  dimension  of  the  parameter  vector. 
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4.3.  LTP-BVP  for  Problem  PI.  We  employ  a backward- 
forward  integration  scheme  in  combination  with  the  method  of 
particular  solutions  (Ref.  7).  The  technique  requires  the 
execution  of  q + 1 independent  sweeps  of  the  differential  sys- 
tem (45) -(52),  each  characterized  by  a different  value  of  the 
multiplier  u . Note  that  Eq.  (46)  is  deleted  and  that  the  simpli- 
fications (7) -(8)  are  invoked  in  the  remaining  equations. 

The  generic  sweep  is  started  by  assigning  particular 
values  to  the  components  of  u;  then,  the  multiplier  A(l)  is 
obtained  from  (50)  . Next,  Eq.  (49)  is  integrated  backward  to 
obtain  the  function  A(t).  With  A (t)  known,  Eq.  (51)  is 
employed  to  obtain  B(t),  and  Eq.  (52)  is  employed  to  compute 
C.  Then,  A ( t)  is  obtained  by  forward  integration  of  (45) 
subject  to  the  initial  condition  (47) . In  this  way,  the 
sweep  is  completed:  for  the  arbitrary  value  assigned  to  p,  it 
leads  to  the  satisfaction  of  all  of  the  equations  of  the 
system  (45)— (52) , except  Eq.  (48). 

In  order  to  satisfy  Eq.  (48)  and  because  the  system 
(45)  — (52)  is  nonhomogeneous , q+1  independent  sweeps  must  be 
executed  employing  q+1  different  multiplier  vectors 
i=l,...,q  + l.  The  first  q sweeps  are  performed  by  choosing 
the  vectors  p^,...,y  to  be  the  columns  of  the  identity 
matrix  of  order  q.  The  last  sweep  is  executed  by  choosing 
to  be  the  null  vector.  As  a result,  one  generates  the 


■nr 


23 


AP-16 


functions  and  multipliers 

Ai(t),  Bi(t),  C.  , Xi(t)>  ui#  i = 1, . . . . ,q + 1.  (55) 

Now,  we  introduce  the  q+1  undetermined,  scalar  constants  k^ 
and  form  the  linear  combinations 

A(t)  = EkiAi(t)  , B(t)  = ZkiBi  (t)  , C=ZkiCi,  (56) 


I 

I 

I 

I 

I 


A(t)  = ZkiAi(t)  , U = ZkiMi  , (57) 

where  the  summations  are  taken  over  the  index  i.  The  q + 1 coef- 
ficients k^  are  obtained  by  forcing  the  linear  combinations  (56) 
to  satisfy  Eq.  (48),  together  with  the  normalization  condition  (Ref.  7) 

Zkt  = 1 . (58) 

Once  the  constants  k^  are  known,  the  solution  of  the  LTP-BVP 
(45)  — (52)  is  given  by  (56)  — (57) . 

4.4.  LTP-BVP  for  Problem  P2.  We  employ  a forward  in- 
tegration scheme  in  combination  with  the  method  of  particular 
solutions  (Ref.  7) . The  technique  requires  the  execution  of 
n + p+1  independent  sweeps  of  the  differential  system  (45M52), 
each  characterized  by  a different  value  of  the  (n  + p) -vector 
a,  whose  components  are  the  n components  of  the  initial  multi- 
plier A ( 0 ) and  the  p components  of  the  parameter  C, 

The  generic  sweep  i3  started  by  assigning  particular 

S ■<  ' /.  I 


values  to  the  components  of  a,  that  is,  the  components  of  the 


%;«*  « <, 
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vectors  A(0)  and  C.  Note  that  A(0)  is  known,  because  of  (47). 
Then,  A ( t)  and  A (t) , together  with  B(t)  and  p(t),  are  obtained 
by  forward  integration  of  (45)  and  (49),  subject  to  (46)  and 
(51).  Note  that,  at  each  time  station  t,  Eqs . (46)  and  (51) 
constitute  a system  of  m + k linear  relations  in  which  the  un- 
knowns are  the  m+k  components  of  the  vectors  B(t)  and  p(t). 
For  this  system  to  have  a unique  solution,  the  following  dis- 

g 

equation  must  hold: 

det  [SuSu  ] * °'  0 < t < 1.  (59) 

As  a result  of  the  procedure,  the  sweep  is  completed:  for  the 
arbitrary  value  assigned  to  a,  it  leads  to  the  satisfaction 
of  all  of  the  equations  of  the  system  (45)-(52),  except  Eqs. 
(48),  (50),  (52). 

In  order  to  satisfy  Eqs.  (48),  (50),  (52)  and  because  the 
system  (45)-(52)  is  nonhomogeneous , n + p+1  independent  sweeps 
must  be  executed  employing  n + p+1  different  vectors 
i = 1 , . . . , n + p + 1 . The  first  n + p sweeps  are  performed  by 


g 

Disequation  (59)  is  obtained  from  (46)  and  (51)  after  elim- 
ination of  B ( t) . The  resulting  linear  equation  in  p(t) 
admits  a unique  solution  providing  (59)  is  satisfied. 
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choosing  the  vectors  o a to  be  the  columns  of  the 

identity  matrix  of  order  n + p.  The  last  sweep  is  executed  by 

choosing  a , . . to  be  the  null  vector.  As  a result,  one 

n + p+1 

generates  the  functions  and  multipliers 

Aj_(t),  Bi(t),  Ci#  A^t),  Oj^Ct),  i = l, ,n  + p + l.  (60) 

Now,  we  introduce  the  n + p+1  undetermined,  scalar  con- 
stants k^  and  form  the  linear  combinations 


A ( t ) = ZkiA;.  (t)  , 

B ( t ) = Zk.B. (t) , 

C = Zk.C.  , 

1 X 

(61) 

1 1 

A(t)  = ZkiAi(t) , 

Pit)  = lkiPi(t)  , 

(62) 

where  the  summations  are  taken  over  the  index  i.  The  n + p+1 
coefficients  k^  and  the  q components  of  the  multiplier  u are 
obtained  by  forcing  the  linear  combinations  (61) -(62)  to 
satisfy  Eqs.  (48),  (50),  (52)  together  with  the  normalization 

condition  (Ref.  7) 

Zki=l.  (63) 


Once  the  constants  k^  are  known,  the  solution  of  the  LTP-BVP 
(45) - (52)  is  given  by  (61)-(62). 

4.5.  Computational  Effort.  Each  sweep  involves 

integrating  2n  differential  equations,  that  is,  the  n linear- 
ized state  equations  (45)  and  the  n multiplier  equations  (49) . 


\ 
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Since  q + 1 sweeps  are  involved  in  Problem  Pi  and  n + p + 1 
sweeps  are  involved  in  Problem  P2,  the  amount  of  computational 
work  performed  per  each  iteration,  gradient  or  restorative, 
is  proportional  to  the  factor: 


Problem  Pi, 

w=2n(q+l); 

(64) 

Problem  P2, 

w = 2n  (n  + p + 1)  . 

(65) 

4.6.  Remark . For  both  the  gradient  phase  and  the  res- 
toration phase,  a linear,  two-point  boundary-value  problem 
must  be  solved.  Once  the  constants  k^  are  known,  the  compo- 
site solution  is  obtained  via  (56) — (57)  or  (61)— (62) . A 
drawback  of  this  procedure  is  that  the  q + 1 particular  solu- 
tions (55)  or  the  n + p + 1 particular  solutions  (60)  must  be 
stored  at  N + 1 time  stations  (here,  N denotes  the  number  of 
integration  subintervals,  so  that  At = 1/N  is  the  magnitude  of 
the  integration  step).  Hence,  a storage  problem  arises  if 
the  system  under  consideration  is  relatively  large,  while  the 
computer  memory  is  relatively  limited. 

This  drawback  can  be  offset  as  follows.  Once  the  con- 
stants k^  are  known,  the  multiplier  u or  the  vector 
a=  [XT(0),CT]T  of  the  composite  solution  is  computed  as 
follows : 
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Problem  Pi, 

M = •••'  kq]  /" 

(66) 

Problem  P2, 

a = [k1,k2,  . . . , kn+p]T- 

(67) 

Then,  a supplementary  sweep  is  executed  according  to  the  procedure 

of  Section  4. 3 or  Section  4.4.  Clearly,  the  total  number  of  sweeps 

increases  by  one,  and  the  computational  work  per  iteration, 

gradient  or  restorative,  becomes  proportional  to  the  factor: 

Problem  Pi, 

w=  2n  (q  + 2)  ; 

(68) 

Problem  P2, 

w=  2n  (n  + p + 2 ) . 

(69) 

In  conclusion,  use  of  this 

supplementary  sweep  increases 

the 

CPU  time;  nevertheless,  depending  on  the  severity  of  the 
storage  problem,  this  course  of  action  might  be  desirable, 
and  sometimes  essential,  with  certain  systems  and  certain 
computers . 

4.7.  Descent  Properties.  The  functions  A(t),  B(t),  C 
solving  Eqs.  (45) -(52)  are  such  that  the  following  first- 

g 

variation  properties  hold: 


gradient  phase. 

6 1 = 5 J = -aQ; 

(70) 

restoration  phase, 

5P = -2aP . 

(71) 

In  Eq.  (70) , Q denotes  the  error  in  the  optimality  conditions 
(16)  at  the  beginning  of  the  gradient  phase.  Because  of 

^This  can  be  3een  by  substitution  of  (45) -(52)  into  (25) -(27). 
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(37)-(40),  the  error  Q reduces  to 

f1 

Q = \ BTBdt  + CTC  . (72) 

J 0 

In  Eq.  (71)  , P denotes  the  constraint  error  (15)  at  the 
beginning  of  the  restoration  phase. 

Note  that  the  first  variation  properties  (70)  and  (71) 
are  consistent  with  the  requirements  (24).  This  being  the 
case,  it  is  possible  in  principle  to  determine  the  gradient 
stepsize  so  that  the  descent  property  (23-1)  or  (23-2)  is 
enforced  in  the  gradient  phase.  It  is  also  possible  to 
determine  the  restoration  stepsize  so  that  the  descent  pro- 
perty (23-3)  is  enforced  in  the  restoration  phase.  For  the 
details,  see  the  following  section. 
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5 . Determination  of  the  Stepsizes 

5.1.  Gradient  Stepsize.  Suppose  that  the  perturbations 
A ( t) , B ( t)  , C solving  the  LTP-BVP  (45)— (52)  for  the  gradient 
phase1^  are  known.  Since  the  nominal  functions  x(t),  u(t),  it 
are  known,  the  one-parameter  family  of  varied  functions  (22) 
can  be  formed.  After  substitution  of  Eqs . (22)  into  Eqs.  (1), 

(9) , (15) , the  following  functions  of  the  stepsize  are 

obtained : 

I = I (a)  , J = J (a)  , P = P (a)  . (73) 

Then,  a one-dimensional  search  scheme  is  applied  to  (73-2) , 
and  a value  of  the  stepsize  a is  selected  for  which  the  fol- 
lowing relations  are  satisfied: 

J (a)  < J (0)  , P(a)<P*,  t (a)  > 0 , (74) 

where  t is  the  final  time  and  P*  is  a preselected  number,  not 
necessarily  small.  Satisfaction  of  Ineq.  (74-1)  is  possible 
because  of  the  descent  property  of  the  gradient  phase.  Ineq. 
(74-2)  is  introduced  to  prevent  excessive  constraint 


iO 


Therefore,  the  constants  kg  and  kr  are  given  by  Eqs,  (53). 
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violation.  And  Ineq.  (74-3)  is  required  for  problems  with 
free  final  time. 

Prior  to  the  satisfaction  of  (74)  , a scanning  process  is 
employed,  leading  to  the  bracketing  of  the  minimum  point  for 
J (a) . This  operation  is  then  followed  by  a Hermitian  cubic 
interpolation  process  (Ref.  34),  which  is  stopped  whenever  the 
following  relation  is  satisfied:'*'1 

^a(a)i-e3  or  I (ot ) ( 0 ) | < e4,  (75) 

subject  to  an  upper  limit  for  the  number  of  search  steps  N . 
Once  a stepsize  a Q has  been  selected  consistently  with  either 
(75)  or  the  prescribed  upper  limit  for  the  number  of  search 
steps,  Ineqs.(74)  must  be  checked.  If  satisfaction  occurs , 
then  the  stepsize  aQ  is  accepted.  If  any  violation  occurs,  then 
the  stepsize  aQ  must  be  bisected  progressively  until  satisfac- 
tion of  (74)  is  finally  achieved. 

5.2.  Remark.  Alternatively,  the  search  for  the  gradient 
stepsize  can  be  performed  by  replacing  the  augmented  function- 
al J with  the  augmented  penalty  functional  W,  defined  by 

W = J + kP , (76) 


11 

The  symbols  e3  and  e4  denote  small,  preselected  numbers. 
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where  k denotes  the  penalty  constant,  to  be  suitably  chosen. 
After  substitution  of  Eqs.  (22)  into  (76),  the  following 
function  of  the  stepsize  is  obtained: 


W = W (a)  . 


(77] 


Then,  the  combination  of  scanning  process/Hermitian  cubic 
interpolation  process  leading  to  satisfaction  of  (75)  is  re- 
placed by  a combination  of  scanning  process/Hermitian  cubic 

interpolation  leading  to  satisfaction  of  the  following  rela- 
12 

tion : 


|Wa(cx)  | < e5 


or 


I Va)/V0)  1 - e6' 


(78) 


subject  to  an  upper  limit  for  the  number  of  search  steps  Ng. 
Once  a stepsize  aQ  has  been  selected  consistently  with  either 
(78)  or  the  prescribed  upper  limit  for  the  number  of  search 
steps,  Ineqs.  (74)  must  be  checked.  If  satisfaction  occurs, 
then  the  stepsize  aQ  is  accepted.  if  any  violation  occurs, 
then  the  stepsize  aQ  must  be  bisected  progressively  until 
satisfaction  of  (74)  is  finally  achieved. 


12 


The  symbols  e5  and  denote  small,  preselected  numbers. 
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5.3.  Restoration  Stepsize.  Suppose  that  the  perturba- 
tions A ( t)  , B ( t ) , C solving  the  LTP-BVP  (45)  — (52)  for  the 
restoration  phase13  are  known.  Since  the  nominal  functions 
x(t),  u(t),  tt  are  known,  the  one-parameter  family  of  varied 
functions  (22)  can  be  formed.  For  this  one-parameter  family, 
the  constraint  error  (15)  becomes  a function  of  the  form 


P = P (a)  . 

Then,  the  stepsize  a must  be  selected  so  that  the  following 
relations  are  satisfied: 


(79: 


P (a)  < P (0)  , 


t (a)  > 0, 


(80) 


Satisfaction  of  Ineq.  (30-1)  is  possible  because  of  the 
descent  property  of  the  restoration  phase.  Ineq.  (80-2)  is 
required  for  problems  with  free  final  time. 

In  order  to  achieve  satisfaction  of  (80) , a bisection 
process  is  applied  to  the  restoration  stepsize  a,  starting 
from  the  reference  stepsize  aQ=l.  This  reference  stepsize 
has  the  property  of  yielding  one-step  restoration  for  the  case 
where  the  constraints  (2) -(5)  are  linear. 


13 


Therefore,  the  constants  k and  kr  are  given  by  Eqs.  (54), 


m 
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5.4.  Iterative  Procedure  for  the  Restoration  Phase. 

The  descent  property  (71)  of  the  restoration  phase  guarantees 
satisfaction  of  Ineq.  (30-1)  at  the  end  of  any  iteration,  but 
not  satisfaction  of  Ineq.  (19-1).  Therefore,  the  restoration 
algorithm  must  be  employed  iteratively  until  Ineq.  (19-1)  is 
satisfied.  At  this  point,  the  restoration  phase  is  termina- 
ted. 


5.5.  Descent  Property  of  a Cycle.  A descent  property 
exists  for  a complete  gradient-restoration  cycle  under  the 
assumption  of  small  stepsizes.  Let  a ^ denote  the  gradient 
stepsize  and  ar  the  restoration  stepsize.  Simple  manipula- 
tions, omitted  for  the  sake  of  brevity,  show  that  th^  gradient 

corrections  are  of  O(oig),  while  the  restoration  corrections 
2 

are  of  0(a  a ).  Hence,  for  a sufficiently  small,  the  restor- 
* a 9 

ation  corrections  are  negligible  with  respect  to  the  gradient 
corrections.  Therefore,  the  restoration  phase  preserves  the 
descent  property  of  the  gradient  phase. 

More  specifically,  let  I,  I,  I denote  the  values  of  the 
functional  (1)  at  the  beginning  of  the  gradient  phase,  at  the 
end  of  the  gradient  phase,  and  at  the  end  of  the  subsequent 
restoration  phase.  Note  that  I and  I are  not  comparable, 
since  the  constraints  are  not  satisfied  to  the  same  accuracy. 
On  the  other  hand,  I and  I are  comparable,  and  the  gradient 
stepsize  can  be  selected  so  that 

A 

I < I . 


(81) 


34 


AP-16 


This  inequality  constitutes  the  descent  property  of  a complete 
gradient-restoration  cycle.  In  order  to  enforce  it,  one  pro- 
ceeds as  follows.  At  the  end  of  the  restoration  phase,  one 
must  verify  Ineq.  (81) . If  it  is  satisfied,  the  next  gradient 
phase  is  started;  otherwise,  the  previous  gradient  stepsize 
is  bisected  as  many  times  as  needed  until,  after  restoration, 
Ineq.  (81)  is  satisfied. 


,SS* 


! 


i 

i 


i 

! 

I 

! 


"H  ✓ 


P 


35 


AP-16 


6 . Summary  of  the  Algorithm 

This  algorithm  includes  cycles  composed  of  a gradient 
phase  and  a restoration  phase.  The  objective  of  each  cycle 
is  to  decrease  the  functional  I so  that  Ineq.  (81)  is  satis- 
fied, while  the  constraints  are  satisfied  to  a predetermined 
accuracy  (19-1) . 

6.1.  Gradient  Phase.  This  phase  involves  a single 
iteration,  and  its  objective  is  to  decrease  the  augmented 
functional  J,  while  the  constraints  are  satisfied  to  first 
order.  The  gradient  phase  can  be  summarized  as  follows. 

(a)  Assume  nominal  functions  x(t),  u(t),rr  which  satisfy 
the  constraints  (2) -(5)  within  the  preselected  accuracy  (19-1). 

(b)  For  the  nominal  functions,  compute  the  vectors  f ^ , 

f , f and  the  matrices  <J>  , <p  , <p  , S , S , S along  the 
interval  of  integration.  At  the  final  point,  compute  the 
vectors  g , g and  the  matrices  . 

(c)  Solve  the  LTP-BVP  (45)- (52) , with  constants  k^  and 
kr  given  by  Eqs.  (53) , using  the  method  of  particular  solu- 
tions. In  this  way,  obtain  the  functions  A(t) , B(t),  C and 
the  multipliers  \ (t) , p (t) , y. 

(d)  Using  the  functions  in  (c) , compute  the  gradient 
stepsize  by  a one-dimensional  search  on  the  augmented  fun- 
ctional J(a)  until  satisfaction  of  Ineq.  (75)  occurs.  Then, 
bisect  the  resulting  stepsize  aQ  (if  necessary) , until 


J 
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satisfaction  of  Ineqs.  (74)  occurs. 

(e)  Once  the  gradient  stepsize  is  known,  compute  the 
varied  functions  x(t),  u(t),  tt  with  Eqs.  (22). 

6.2.  Restoration  Phase.  This  phase  involves  one  or 
more  iterations,  and  its  objective  is  to  reduce  the  constraint 
error  P,  until  satisfaction  of  (19-1)  occurs.  Within  a single 
iteration,  the  objective  is  to  decrease  the  constraint  error 
to  a level  compatible  with  Ineq.  (23-3) , while  the  norm 
squared  of  the  variations  of  the  control  and  the  parameter 
is  minimized. 

The  nominal  functions  x(t),  u(t),  tt  are  chosen  as  fol- 
lows: for  the  first  restorative  iteration,  the  nominal  fun- 

ctions are  identical  with  the  varied  functions  obtained  at 
the  end  of  the  previous  gradient  iteration;  for  any  subsequent 
restorative  iteration,  the  nominal  functions  are  identical 
with  the  varied  functions  obtained  at  the  end  of  the  previous 
restorative  iteration.  With  this  understanding,  the  resto- 
ation  phase  can  be  summarized  as  follows. 

(a)  Assume  nominal  functions  x(t),  u(t),  n which  satisfy 
condition  (3),  but  violate  at  least  one  of  conditions  (2)  and 
(4)  - (5)  . 

(b)  For  the  nominal  functions,  compute  the  vectors 

(x  - <J>)  , S and  the  matrices  <j>  , <J>u,  Sx,  Su,  along  the 

interval  of  integration.  At  the  final  point,  compute 
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the  vector  p and  the  matrices  p , p^. 

(c)  Solve  the  LTP-BVP  ( 4 5 ) - ( 52 ),  with  constants  k and 

g 

k^  given  by  Eqs.  (54),  using  the  method  of  particular  solu- 
tions. In  this  way,  obtain  the  functions  A(t),  B(t),  C and 
the  multipliers  X (t) , p(t),  y • 

(a)  Using  the  functions  in  (c) , compute  the  restoration 
stepsize  by  a one-dimensional  search  on  the  constraint  error 
P(a).  To  this  effect,  perform  a bisection  process  on  a, 
starting  from  aQ=l,  until  Ineqs.  (80)  are  satisfied. 

(e)  Once  the  restoration  stepsize  is  known,  compute  the 
varied  functions  x(t),  u(t),  tt  with  Eqs.  (22). 

(f)  Verify  whether  the  varied  functions  in  (e)  satisfy 
Ineq.  (19-1) . If  this  is  the  case,  the  restoration  phase  is 
terminated.  Otherwise,  return  to  (a)  and  continue  the  pro- 
cess until  satisfaction  of  (19-1)  occurs. 

6.3.  Gradient-Restoration  Cycle.  After  the  restoration 
phase  is  completed,  verify  whether  Ineq.  (81)  is  satisfied. 

If  this  is  the  case,  start  the  next  cycle  of  the  sequential 
gradient-restoration  algorithm.  If  not,  return  to  the  previous 
gradient  phase  and  reduce  the  gradient  stepsize  (using  a bisection 
process)  until,  after  restoration,  Ineq.  (81)  is  satisfied. 

6.4.  Computational  Considerations.  Here,  special 
conditions  relevant  to  the  computer  implementation  of  the 
sequential  gradient-restoration  algorithm  are  presented. 
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Starting  Condition.  The  present  algorithm  can  be 
started  with  nominal  functions  x(t),  u(t),  tt  satisfying  con- 
dition (3)  and  violating  none,  some,  or  all  of  conditions  (2) 
and  (4) -(5).  If  the  nominal  functions  are  such  that  Ineq. 
(19-1)  is  violated,  the  algorithm  starts  with  a restoration 
phase;  hence,  the  first  cycle  is  a half  cycle,  involving  a 
restoration  phase  only.  On  the  other  hand,  if  the  nominal 
functions  are  such  that  Ineq.  (19-1)  is  satisfied,  the  algor- 
ithm starts  with  a gradient  phase;  hence,  the  first  cycle  is 
a complete  gradient-restoration  cycle. 

Bypassing  Condition.  At  the  end  of  the  gradient  phase 
of  any  cycle,  the  constraint  error  P must  be  computed.  If 
Ineq.  (19-1)  is  violated,  a restoration  phase  is  started. 
Otherwise,  the  restoration  phase  is  bypassed,  and  the  next 
gradient  phase  of  the  algorithm  is  started. 

Stopping  Conditions.  For  the  restoration  phase  taken 
individually,  convergence  is  achieved  whenever  Ineq.  (19-1)  is 
satisfied.  For  the  sequential  gradient-restoration  algorithm 
taken  as  a whole,  convergence  is  achieved  whenever  Ineqs. 

(19-1)  and  (19-2)  are  satisfied  simultaneously. 
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7 . Experimental  Conditions 

In  order  to  evaluate  the  theory,  several  examples  were 
solved.  The  sequential  gradient-restoration  algorithms 
associated  with  Problems  Pi  and  P2  were  programmed  in 
FORTRAN  IV,  and  the  numerical  results  were  obtained  in  double- 
precision arithmetic. 

Computations  were  performed  at  Rice  University  using  an 
IBM  370/155  computer.  For  each  example,  the  interval  of  in- 
tegration was  divided  into  100  steps.  The  differential 
equations  were  integrated  using  Hamming's  modified  predictor- 
corrector  method  with  a special  Runge-Kutta  starting  procedure 
(Ref.  35).  The  definite  integrals  I,  J,  P,  Q were  computed 
using  a modified  Simpson's  rule.  The  method  of  particular 
solutions  (Ref.  7)  was  used  to  solve  the  linear,  two-point 
boundary-value  problems  associated  with  both  the  gradient 
phase  and  the  restoration  phase. 

7.1.  Convergence  Conditions.  The  parameters  e^, 

' ~ 14 

e4  appearing  in  Ineqs.  (19)  and  (75)  were  set  at  the  levels 

£.  = E-08 , e.  = E-04,  e . = E-03  . (82) 

1 2 4 


;j  

I^The  symbol  E + ab  stands  for  lO-3*3 

!;  ! 
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The  tolerance  level  (82-1)  characterizes  the  restoration 
phase;  the  tolerance  levels  (82-1)  and  (82-2),  employed  in 
combination,  characterize  the  algorithm  as  a whole;  and  the 
tolerance  level  (82-3)  characterizes  the  one-dimensional 
search  for  the  gradient  stepsize. 

7.2.  Safeguards . For  the  gradient  phase,  the  parameter 
P*  appearing  in  Ineq.  (74-2)  was  set  at  the  level 

P*  = 10  . (83) 

The  tolerance  level  (83)  limits  the  constraint  violation 
which  is  permissible  during  the  gradient  phase.  Also  for  the 
gradient  phase,  the  number  of  Hermitian  search  steps  required 
to  satisfy  Ineq.  (75)  was  subject  to  the  upper  bound 


Ng  < 5 . (34) 

7.3.  Nonconvergence  Conditions.  The  sequential 
gradient-restoration  algorithms  were  programmed  to  stop  when- 
ever violation  of  any  of  the  following  inequalities  occurred:^ 


I 


I 

I 

I 


15 


Inequality  (87)  is  characteristic  of  the  IBM  370/155  computer. 
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' < 30,  N < 100, 

N < 10, 
r - ' 

(85) 

.<10,  N.  < 10, 

bg  - ' br  - 

Nbc±5' 

(86) 

M < 0.83  E + 75  . 

(87) 

Here,  Nc  is  the  number  of  cycles,  N is  the  total  number  of 
iterations,  Nr  is  the  number  of  restorative  iterations  per 
cycle,  is  the  number  of  bisections  of  the  gradient  step- 

size  required  to  satisfy  Ineqs.  (74),  is  the  number  of 

bisections  of  the  restoration  stepsize  required  to  satisfy 
Ineqs.  (80) , Nfac  is  the  number  of  bisections  of  the  gradient 
stepsize  required  to  satisfy  Ineq.  (81) , and  M is  the  modulus 
of  any  of  the  quantities  employed  in  the  algorithm. 
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8 . Numerical  Examples,  Problem  PI 

In  this  section,  four  numerical  examples  are  described 
employing  scalar  notation.  In  particular,  the  symbols 
x^(t),  i = l,...,n,  denote  the  components  of  the  state;  the 
symbols  u^(t),  i = l,...,m,  denote  the  components  of  the  con- 
trol; and  the  symbols  tt^,  i=l,...,p,  denote  the  components 
of  the  parameter. 

For  all  of  the  examples,  a time  normalization  is  used  in 
order  to  simplify  the  numerical  computations.  Specifically, 
the  actual  time  0 is  replaced  by  the  normalized  time 

t-e/T,  (88) 

which  is  defined  in  such  a way  that  t = 0 at  the  initial  point 
and  t = 1 at  the  final  point.  The  actual  final  time  t,  if  it 
is  free,  is  regarded  as  a component  of  the  vector  parameter  tt 
to  be  optimized.  In  this  way,  an  optimal  control  problem  with 
variable  final  time  is  converted  into  an  optimal  control  prob- 
lem with  fixed  final  time. 

Concerning  the  convergence  history,  the  terminology  is 

as  follows:  N denotes  the  cycle  number,  N is  the  number  of 
c g 

gradient  iterations  per  cycle,  Nr  is  the  number  of  restorative 
iterations  per  cycle,  N is  the  total  number  of  iterations, 

P is  the  constraint  error,  Q is  the  error  in  the  optimality 
conditions,  and  I the  value  of  the  functional  being  minimized. 
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Example  8.1.  This  example  involves  (i)  a quadratic 
functional,  (ii)  nonlinear  differential  equations,  (iii) 
boundary  conditions  of  the  fixed  endpoint  type,  and  (iv) 
fixed  final  time  x = l: 


f1  . . 

I = \ (1  + xj  + x2  + uf)dt  , 

(39) 

Jo 

2 

= u^  - x2  , x2  = u-j^  - x1x2  , 

(90) 

xx(0)  = 0,  x2(0)  = 1 , 

(91) 

x1  (1)  =1,  x2  (1)  = 2 . 

(92) 

The  assumed  nominal  functions  are: 

x1(t)=t,  x2(t)=l  + t,  u1(t)=l.  (93) 

The  numerical  results  are  given  in  Tables  1-2.  Convergence  to 
the  desired  stopping  condition  occurs  in  Nc  = 3 cycles  and 
N=7  iterations,  which  include  2 gradient  iterations  and  5 
restorative  iterations. 

Example  8.2.  This  example  involves  (i)  a nonquadratic 
functional,  (ii)  nonlinear  differential  equations,  (iii) 
boundary  conditions  of  the  fixed  endpoint  type,  and  (iv)  fixed 
final  time  t*  1: 
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r1 

I=\  (-2cos  u.  ) dt  , 

(94) 

Jo 

x^  = 2sin  u^  - 1 , 

x2  = xx  , 

(95) 

xx  (0)  =0  , 

x2(0)  = 0, 

(96) 

xx  (1)  =0  , 

x2(l)  = 0.3  . 

(97) 

The  assumed  nominal  functions  are: 

x1(t)=0,  x2(t)«0.3t,  ux(t)=0.  (98) 

The  numerical  results  are  given  in  Tables  3-4.  Convergence 
to  the  desired  stopping  condition  occurs  in  Nc  = 6 cycles  and 
N = 13  iterations,  which  include  5 gradient  iterations  and  8 
restorative  iterations. 

Example  8.3.  This  example  is  a minimum  time  problem  and 
involves  (i)  a linear  functional,  (ii)  nonlinear  differential 
equations,  (iii)  boundary  conditions  of  the  fixed  final  state 
type,  and  (iv)  free  final  time  r.  After  setting  = the 
problem  is  as  follows: 


I = 


(99) 


X1  = 1T1U1' 


. . 2 2, 
x2  = 'rr1  - x^)  , 


(100) 


' v 


ki'T  '!• 
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xx(0)  = 0, 

X 

to 

0 

II 

0 

(101) 

xL (1)  = 1 , 

0 

II 

r— I 

CM 

X 

(102) 

The  assumed  nominal  functions 

are : 

x^  (t)  = t,  x2  (t)  = 0 , 

u1 (t)  =1, 

tt^  = 1 . 

(103) 

The  numerical  results  are  given  in  Tables  5-6.  Convergence 
to  the  desired  stopping  condition  occurs  in  Nc  = 3 cycles  and 
N = 7 iterations,  which  include  2 gradient  iterations  and  5 
restorative  iterations. 

Example  3.4.  This  example  is  a minimum  time  problem  and 
involves  (i)  a linear  functional,  (ii)  nonlinear  differential 
equations,  (iii)  components  of  the  final  state  partly  given 
and  partly  free,  and  (iv)  free  final  time  t.  After  setting 
= the  problem  is  as  follows: 

I = tt1,  (104) 

xL  = tt^x^cos  ux  , *2~  irix3s^n  > X3  = rr^sin  u^  , (105) 

Xj^OJ-O,  x2(0)=0,  x3(0)=0,  (106) 

xx(l)  = 1 . (107) 
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The  assumed  nominal  functions  are: 

x1(t)=t,  x2(t)=0,  x3(t)  = 0,  u^(t)=l,  71^=1.  (103) 

The  numerical  results  are  given  in  Tables  7-8-  Convergence  to 
the  desired  stopping  condition  occurs  in  Nc  = 4 cycles  and 
N = 12  iterations,  which  include  3 gradient  iterations  and  9 
restorative  iterations. 
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9 . Numerical  Examples,  Problem  P2 

In  this  section,  four  numerical  examples  are  described 
employing  scalar  notation.  The  symbols  used  for  Problem  P2 
are  the  same  as  those  employed  for  Problem  Pi.  The  time 
normalization  (88)  is  employed. 

Example  9.1.  This  example  involves  (i)  a quadratic 
functional,  (ii)  a nonlinear  differential  equation,  (iii)  a 
state  inequality  constraint  of  the  first  order^,  (iv)  bound- 
ary conditions  of  the  fixed  endpoint  type,  and  (v)  fixed  final 
time  T = 1 : 


1 = \ (xj  + uj)dt  , 


X1  = X1_U1  ' 


(110) 


xx  - 0.9  > 0 , 


Xjl  (0)  =1  , 


xx(l)  = 1 . 


(113) 


This  means  that  the  first  time  derivative  of  the  left  hand 
side  of  Ineq.  (Ill)  contains  the  control  explicitly. 
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Upon  introducing  the  auxiliary  state  variable  x2  and  the 
auxiliary  control  variable  U2  defined  by  (Ref.  25) 


x 


1 


0.9 


(114) 


we  replace  the  inequality  constrained  problem  (109) -(113)  with 
the  following  equality  constrained  problem: 

I = ( (xj*  + u*)dt  , (115) 


JO 

• 2 • 

X1  = x£  - u1  , x2  = u2  , (116) 

x^  - U1  - 2x2u2  = 0 , (117) 

xx(0)  = 1 , x2(0)  = /(0.1)  , (118) 

x1(l)  = 1 . (119) 


The  assumed  nominal  functions  are: 

x1(t)=l,  x2  (t)= /(0.1)  , u1(t)=l,  u2(t)=l.  (120) 

The  numerical  results  are  given  in  Tables  9-10.  Convergence 
to  the  desired  stopping  condition  occurs  in  Nc  = 5 cycles  and 
N = 12  iterations,  which  include  4 gradient  iterations  and  8 
restorative  iterations. 
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Example  9.2.  This  example  involves  (i)  a quadratic 
functional,  (ii)  linear  differential  equations,  (iii)  a state 
inequality  constraint  of  the  second  order17,  (iv)  boundary 
conditions  of  the  fixed  endpoint  type,  and  (v)  fixed  final 


I = 

(121) 

Xx  — ^2  f 

Jo 

x2  = ui  , 

(122) 

0.15 

0 

A 1 

H 

X 

1 

(123) 

xx  (0)  = 0 , 

x2(0)  = 1 , 

. (124) 

xl  (1)  =0  , 

x2(l)  = -1  . 

(125) 

Upon  introducing  the  auxiliary  state  variables  x3 , x4  and  the 
auxiliary  control  variable  u^  defined  by  (Ref.  25) 

2 

0.15-x1*x3,  x3  = x4  , x4  = u2  , (126) 


This  means  that  the  second  time  derivative  of  the  left-hand 
side  of  Ineq.  (123)  contains  the  control  explicitly,  while 
this  is  not  the  case  with  the  first  time  derivative. 
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we  replace  the  inequality  constrained  problem  (121) -(125)  with 
the  following  equality  constrained  problem: 

I = ( uj  dt  , (127) 


Xx  — X2  / 


x2  " U1  ' 


x3  = x4  , 


x4  - u2  / 


(123) 


u,  + 2x,u„  + 2x“T  = 0 , 
1 3 2 4 


XL(1)  = 0, 


x2  (1)  = -1 


The  assumed  nominal  functions  are: 


x1(t)=0,  x2(t)=l-2t. 


x4  (t)  = (2t  - l)//(0.60)  , u1(t)=l,  u2(t)=0 


(129) 


x1(0)=0,  x2  (0)  = 1 , x3(0)=  /(0.15)  , x4  (0)  = -i//(0. 60)  , (130) 


(131) 


x- (t)  = / (0.15)  (1  - 2t) , (132) 


(133) 


The  numerical  results  are  given  in  Tables  11-12.  Convergence 
to  the  desired  stopping  condition  occurs  in  Nc  = 8 cycles  and 
N = 16  iterations,  which  include  7 gradient  iterations  and  9 
restorative  iterations. 

Example  9.3.  This  example  is  a minimum  time  problem  and 
involves  (i)  a linear  functional,  (ii)  nonlinear  differential 
equations,  (iii)  a state-derivative  inequality  constraint, 


- . ; 
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(iv)  boundary  conditions  of  the  fixed  final  state  type,  and  (v) 
free  final  time  x.  After  setting  the  problem  is  as 

follows : 


I — Tf  ^ , 


(134) 


X1  = *1*1  ' 


• 2 2 

X2  = U1  ^U1  ~ xi  ~ 0 • 5)  > 


(135) 


x2/ir  + 0 . 5 > 0 , 


(136) 


x1  (0)  = 0 , x2  (0)  = 0 , 


(137) 


x^l)  = 1 , 


x2  (1)  = -n/4  . 


(138) 


Upon  introducing  the  auxiliary  control  variable  u2  defined  by 


(Ref.  25) 


x A + 0 . 5 - u2  = 0 , 


(139) 


we  replace  the  inequality  constrained  problem  (134) -(138)  with 
the  following  equality  constrained  problem: 


I = xr  x , 


(140) 


X1  = 11  iui  ' 


2 2 

x2  ■ (u^  - xl  - 0.5)  , 


(141) 


2 2 2 . 
u1  - x1  - u2  = 0 , 


(142) 


xx(0)  = 0 , 


x2(0)  = 0 , 


(143) 


52 


AP-16 


I 


x1(l)  = 1 , 


x2  (1)  = - tt/4  . 


(144) 


The  assumed  nominal  functions  are: 

x1(t)=t,  x2  (t)  =- ( tt/4)  t,  u1(t)=l,  u2(t)=l,  17^=1.  (145) 

The  numerical  results  are  given  in  Tables  13-14.  Convergence 
to  the  desired  stopping  condition  occurs  in  Nq  = 6 cycles  and 
N = 14  iterations,  which  include  5 gradient  iterations  and  9 
restorative  iterations. 

Example  9.4.  This  example  involves  (i)  a quadratic 
functional,  (ii)  linear  differential  equations,  (iii)  a control 
inequality  constraint,  (iv)  boundary  conditions  of  the  fixed 
endpoint  type,  and  (v)  fixed  final  time  t: 

*1 


1=1  (1  + x^  + x2  + u^) dt  , 

Jo 

(146) 

= ul  - x2  , 

*2  = ui  “ , 

(147) 

6 - u^  > 0 , 

(148) 

(0)  = 0 > 

x2l0)  = 1 , 

(149) 

5^(1)  = 1 , 

x2(.l)  = 2 . 

(150) 

Upon  introducing  the  auxiliary  control  variable  u2  defined  by 


I 

t 
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S 


6 - u^  - 

•Uj-O  , 

(151) 

we  replace  the  inequality  constrained  problem  (146) -(150) 

with  the  following  equality  constrained  problem: 

f1 

1=  '1  + xl 

+ + u^) dt  , 

(152) 

J 0 

*1  = ui  * x2  ' 

x2  = ui  - 2x1  , 

(153) 

6 - u^  - 

• u2  = 0 , 

(154) 

x1(0)  = 0 , 

x2  (0)  =1  , 

(155) 

xx(l)  = 1 , 

x2  (1)  = 2 . 

(156) 

The  assumed  nominal  functions  « 

are : 

x^ (t)  = 5t  - 4t^  , 

x2  ( t)  = 1+  5t  - 4 1^  , 

(157) 

uL(t)  = 6 (1  - t)  , 

u2 (t)  = 2t  . 

(158) 

The  numerical  results  are  given  in  Tables  15-16,  Convergence 
to  the  desired  stopping  condition  occurs  in  Nc  = 11  cycles  and 
N =*  24  iterations,  which  include  10  gradient  iterations  and  14 
restorative  iterations. 
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Table  1.  Convergence  history,  Example  8.1. 


1 

I 

I 


I 

I 

t 

I 

| 

I 

I 


N 

c 

N 

g 

N 

r 

N 

P 

Q 

I 

0 

0 

0 

0 

0.72E+01 

1 

0 

4 

4 

0.32E-10 

0. 97E+00 

33.67701 

2 

l 

1 

6 

0.84E-13 

0.50E-02 

33.46606 

3 

l 

0 

7 

0.51E-09 

0.41E-04 

33.46484 

Table  2. 

Converged 

solution.  Example 

CO 

t 

X1 

x2 

U1 

0.0 

0.0000 

1.0000 

-8.3428 

0.1 

-0.7862 

0.2773 

-6.3676 

0.2 

-1.3011 

-0.2366 

-3.8632 

0.3 

-1.5837 

-0.5625 

-1.4845 

0.4 

-1.6735 

-0.7169 

0.4682 

0.5 

-1.6003 

-0.7107 

1.9931 

0.6 

-1.3780 

-0.5437 

3.2522 

0.7 

-1.0080 

-0.2055 

4.4920 

0.8 

-0.4877 

0.3179 

6.0526 

0.9 

0.1807 

1.0416 

8.4996 

1.0 

1.0000 

2.0000 

13.0496 

T = 

1.00000 

yJt.’V-v.  \ 
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Table  3.  Convergence  history.  Example  8.2. 


N 

c 

N 

g 

N 

r 

N 

P 

Q 

I 

0 

0 

0 

0 

0.10E+01 

1 

0 

4 

4 

0.17E-08 

0.67E+00 

-1.11665 

2 

1 

2 

7 

0.17E-11 

0.34E-01 

-1.16519 

3 

1 

1 

9 

0.11E-10 

0.30E-02 

-1.16923 

4 

1 

1 

11 

0.58E-15 

0.64E-03 

-1.16950 

5 

1 

0 

12 

0.20E-08 

0.18E-03 

-1.16961 

6 

1 

0 

13 

0.36E-08 

0.50E-04 

-1.16964 

Table  4 . 

Converged 

solution,  Example 

8.2. 

t 

X1 

x2 

U1 

0.0 

0.0000 

0.0000 

1.3333 

0.1 

0.0937 

0.0047 

1.3049 

0.2 

0.1856 

0.0186 

1.2609 

0.3 

0.2742 

0.0417 

1.2005 

0.4 

0.3575 

0.0733 

1.1131 

0.5 

0.4309 

0.1128 

0.9784 

0.6 

0.4842 

0.1589 

0.7517 

0.7 

0.4921 

0.2082 

0.3661 

0.8 

0.4141 

0.2544 

-0.1521 

0.9 

0.2381 

0.2877 

-0.6087 

1.0 

0.0000 

0.3000 

-0.8959 

T = 1.00000 
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Table  5.  Convergenge  history,  Example  3.3. 


i 


I 


I 

I 

I 

I 

I 

I 


Nc 

N 

g 

Nr 

N 

P 

Q 

I 

0 

0 

0 

0 

0. 53E+00 

1 

0 

4 

4 

0.74E-16 

0.53E-01 

1.58101 

2 

l 

1 

6 

0.33E-08 

0.13E-03 

1.57080 

3 

l 

0 

7 

0. 28E-08 

0.16E-05 

1.57075 

Table  6. 

Converged 

solution,  Example 

8.3. 

t 

X1 

X2 

U1 

0.0 

0.0000 

0.0000 

0.9997 

0.1 

0.1564 

0.1544 

0.9874 

0.2 

0.3089 

0.2937 

0.9508 

0.3 

0.4538 

0.4043 

0.8907 

0.4 

0.5876 

0.4752 

0.8087 

0.5 

0.7069 

0.4997 

0.7067 

0.6 

0.8087 

0.4752 

0.5875 

0.7 

0.8907 

0.4042 

0.4538 

0.8 

0.9507 

0.2937 

0.3092 

0.9 

0.9874 

0.1544 

0.1572 

1.0 

1.0000 

0.0000 

0.0017 

T = TT1 

» 1.57075 
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I 

I 

I 

I 


I 

I 

I 

I 

I 

I 

I 

I 

! 

I 


Table  7.  Convergence  history,  Example  8.4. 


N 

c 

N 

g 

N 

r 

N 

P 

Q 

I 

0 

0 

0 

0 

0.17E+01 

1 

0 

5 

5 

0.44E-16 

0.25E+00 

1.83370 

2 

1 

2 

8 

0.29E-09 

0.42E-01 

1.78266 

3 

1 

1 

10 

0.39E-09 

0.62E-03 

1.77262 

4 

n 

1 

12 

0.10E-16 

0.92E-05 

1.77245 

Table  8.  Converged  solution,  Example  8.4. 

t 

X1  X2  X3 

U1 

0.0 

0.0000 

0.0000 

0.0000 

1.5708 

0.1 

0.0016 

0.0155 

0.1765 

1.4133 

0.2 

0.0129 

0.0607 

0.3486 

1.2558 

0.3 

0.0425 

0.1311 

0.5121 

1.0984 

0.4 

0.0974 

0.2198 

0.6630 

0.9412 

0.5 

0.1819 

0.3180 

0.7975 

0.7841 

0.6 

0.2976 

0.4161 

0.9123 

0.6270 

0.7 

0.4428 

0.5046 

1.0046 

0.4699 

0.8 

0.6132 

0.5748 

1.0722 

0.3126 

0.9 

0.8018 

0.6196 

1.1132 

0.1549 

1.0 

1.0000 

0.6346 

1.1266 

-0.0033 

x - tv1  - 1.77245 

« 


• -a- 
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Table  9.  Convergence  history,  Example  9.1. 


N 

c 

N 

g 

N 

r 

N 

P 

Q 

I 

0 

0 

0 

0 

0.14E+01 

1 

0 

3 

3 

0. 52E-09 

0.35E+00 

1.83569 

2 

1 

2 

6 

0.15E-16 

0.14E-01 

1.66599 

3 

1 

1 

8 

0.10E-09 

0.24E-03 

1.65742 

4 

1 

1 

10 

0.60E-17 

0.15E-03 

1.65697 

5 

1 

1 

12 

0.96E-18 

0.98E-04 

1.65678 

1 

1 

1 

Table  10.  Converged  solution 

, Example 

9.1. 

1 

t 

| 

X1 

X2 

U1 

U2 

1 

0.0 

1.0000 

0.3162 

1.7482 

-1.1831 

• 0.1 

0.9410 

0.2025 

1.3353 

-1.1104 

0.2 

0.9095 

0.0978 

1.0097 

-0.9324 

0.3 

0.9006 

0.0246 

0.8366 

-0.5177 

, 0.4 

0.9000 

-0.0090 

0.8067 

-0.1865 

0.5 

0.9003 

-0.0177 

0.8104 

-0.0018 

0.6 

0.9000 

-0.0094 

0.8135 

0.1816 

0.7 

0.9005 

0.0238 

0.7864 

0.5158 

0.8 

0.9094 

0.0972 

0.6442 

0.9398 

0.9 

0.9409 

0.2024 

0.4360 

1.1097 

" 1.0 

1.0000 

0.3162 

0.2470 

1.1904 

= 1.00000 
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0 

0 

0 

0 

0. 22E+02 

1 

0 

5 

5 

0.44E-13 

0.11E+00 

6.03009 

2 

1 

1 

7 

0.15E-14 

0.79E-02 

5.93793 

3 

1 

1 

9 

0. 28E-17 

0.20E-02 

5.93016 

4 

1 

1 

11 

0.12E-18 

0.74E-03 

5.92817 

5 

1 

1 

13 

O.J  5E-20 

0.37E-03 

5.92738 

6 

1 

0 

14 

0. 83E-08 

0.20E-03 

5.92687 

7 

1 

0 

15 

0.62E  -^8 

0.12E-03 

5.92661 

8 

1 

0 

16 

0.74E  1 

0.52E-04 

5.92650 

I 

I 

Table  12.  Converged  solution,  Example  9.2. 


t 

X1 

X2 

X3 

X4 

U1 

U2 

0.0 

0.0000 

1.0000 

0.3872 

-1.2909 

-4.4535 

1.4461 

0.1 

0.0793 

0.6045 

0.2657 

-1.1375 

-3.4592 

1.6392 

0.2 

0.1242 

0.3078 

0.1606 

-0.9583 

-2.4746 

1.9862 

0.3 

0.1442 

0.1105 

0.0756 

-0.7301 

-1.4669 

2.6475 

0.4 

0.1496 

0.0145 

0.0175 

-0.4152 

-0.4730 

3.6484 

0.5 

0.1499 

-0.0001 

-0.0045 

-0.0174 

0.0375 

4.1943 

0.6 

0.1497 

-0.0118 

0.0147 

0.4010 

-0.4405 

4.0278 

0.7 

0.1446 

-0.1097 

0.0733 

0.7483 

-1.5262 

2.7693 

0.8 

0.1244 

-0.3098 

0.1599 

0.9688 

-2.4681 

1.8476 

0.9 

0.0794 

-0.6057 

0.2655 

1.1403 

-3.4508 

1.6007 

1.0 

0.0000 

-1.0000 

0.3872 

1.2909 

-4.4354 

1.4228 

t = 1.00000 


Table  13.  Convergence  history,  Example  9.3. 


N 

c 

N 

g 

N 

r 

N 

P 

Q 

I 

0 

0 

0 

0 

0.11E+01 

1 

0 

5 

5 

0.22E-14 

0.21E-01 

1.82848 

2 

l 

2 

8 

0.47E-15 

0. 20E-02 

1.82290 

3 

l 

1 

10 

0.33E-12 

0.55E-03 

1.82245 

4 

l 

1 

12 

0.18E-13 

0.22E-03 

1.82234 

5 

l 

0 

13 

0.60E-08 

0.10E-03 

1.82224 

6 

l 

0 

14 

0.77E-08 

0. 39E-04 

1.82222 

I 


Table  14.  < 

Converged  solution 

, Example 

9.3. 

t 

X1 

X2 

U1 

U2 

0.0 

0.0000 

0.0000 

0.4999 

0.4999 

0.1 

0.0905 

-0.0465 

0.4916 

0.4832 

0.2 

0.1781 

-0.0989 

0.4670 

0.4317 

0.3 

0.2598 

-0.1623 

0.4271 

0.3389 

0.4 

0.3331 

-0.2401 

0.3768 

0.1762 

0.5 

0.4020 

-0.3298 

0.4020 

0.0092 

0.6 

0.4824 

-0.4209 

0.4824 

0.0000 

0.7 

0.5788 

-0.5120 

0.5788 

0.0001 

0.8 

0.6945 

-0.6031 

0.6945 

0.0000 

0.9 

0.8334 

-0.6942 

0.8333 

-0.0007 

1.0 

1.0000 

-0.7853 

0.9996 

i 

0.0008 

T - irL  ■ 1.82222 
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Table  15.  Convergence  history,  Example  9.4. 


N 

c 

N 

g 

N 

r 

N 

P 

Q 

I 

0 

0 

0 

0 

0.38E+01 

1 

0 

4 

4 

0.18E-12 

0. 36E+00 

20.27422 

2 

1 

2 

7 

0. 36E-10 

0.37E-01 

20.19329 

3 

1 

1 

9 

0.79E-10 

0.96E-02 

20.18932 

4 

1 

1 

11 

0.12E-11 

0 . 49E-02 

20.18813 

5 

1 

1 

13 

0.32E-13 

0.20E-02 

20.18760 

6 

1 

1 

15 

0. 34E-14 

0.12E-02 

20.18733 

7 

1 

1 

17 

0. 20E-15 

0.63E-03 

20.18718 

8 

1 

1 

19 

0. 96E-16 

0.50E-03 

20.18707 

9 

1 

1 

21 

0.45E-17 

0.24E-03 

20.18700 

10 

1 

1 

23 

0. 16E-16 

0.34E-03 

20.18693 

11 

1 

0 

24 

0. 29E-08 

0. 71E-04 

20.18688 

I 


Table  16. 

Converged  solution, 

Example  9.4. 

t 

X1 

X2 

U1  U2 

0.0 

0.0000 

1.0000 

6.0000 

0.0000 

0.1 

0.4716 

1.5519 

6.0000 

0.0000 

0.2 

0.8742 

1.9967 

5.3504 

0.8059 

0.3 

1.1410 

2.2746 

4.3155 

1.2978 

0.4 

1.2930 

2.4172 

3.4641 

1.5924 

0.5 

1.3586 

2.4610 

2.7629 

1.7991 

0.6 

1.3598 

2.4346 

2.1838 

1.9535 

0.7 

1.3133 

2.3603 

1.7031 

2.0728 

0.8 

1.2320 

2.2549 

1.3012 

2.1676 

0.9 

1.1253 

2.1315 

0.9618 

2.2445 

1.0 

1.0000 

2.0000 

0.6710 

2.3084 

T = 1.00000 
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1C . Discussion  and  Conclusions 

In  this  paper,  two  members  of  the  family  of  sequential 
gradient-restoration  algorithms  for  the  solution  of  optimal 
control  problems  are  presented.  These  algorithms  are  of  the 
ordinary-gradient  type.  One  is  associated  with  the  solution 
of  Problem  Pi,  Eqs . (l)-(4),  and  the  other  is  associated  with 

the  solution  of  Problem  P2,  Eqs.  ( 1 ) — ( 5 ) . 

Problem  Pi  consists  of  minimizing  a functional  I which 
depends  on  the  n-vector  state  x(t),  the  m-vector  control  u(t), 
and  the  p-vector  parameter  it.  The  state  is  given  at  the  ini- 
tial point.  At  the  final  point,  the  state  and  the  parameter 
are  required  to  satisfy  q scalar  relations.  Along  the  inteval 
of  integration,  the  state,  the  control,  and  the  parameter  are 
required  to  satisfy  n scalar  differential  equations.  Problem 
P2  differs  from  Problem  Pi  in  that  the  state,  the  control, 
and  the  parameter  are  required  to  satisfy  k additional  scalar 
relations  along  the  interval  of  integration. 

The  importance  of  Problems  pi  and  P2  lies  in  the  fact 
that  a large  number  of  problems  of  optimal  control  are  covered 
by  these  formulations  (Refs.  7-34).  In  particular,  Problem 
P2  enlarges  dramatically  the  number  and  variety  of  problems 
of  optimal  control  which  can  be  treated  by  gradient-restoration 
algorithms.  Indeed,  by  suitable  transformations,  almost  every 
known  problem  of  optimal  control  can  be  brought  into  the 


63 


AP-16 


I 

I 

| I 
I 
I 

I 
I 
I 

I 


scheme  of  Problem  P2.  This  statement  applies,  for  instance, 
to  the  following  situations:  (i)  problems  with  control 
equality  constraints,  (ii)  problems  with  state  equality  con- 
straints, (iii)  problems  with  state-derivative  equality  con- 
straints, (iv)  problems  with  control  inequality  constraints, 

(v)  problems  with  state  inequality  constraints,  and  (vi)  prob- 
lems with  state-derivative  inequality  constraints.  For  an 
illustration  of  the  scope  and  range  of  applicability  of 
Problem  P2,  the  reader  is  referred  to  Ref.  19  and  Refs.  25-29. 

The  algorithms  presented  here  include  a sequence  of  two- 
phase  cycles,  composed  of  a gradient  phase  and  a restoration 
phase.  The  gradient  phase  involves  one  iteration  and  is  de- 
signed to  decrease  the  value  of  the  functional  I,  while  the 
constraints  are  satisfied  to  first  order.  The  restoration 
phase  involves  one  or  more  iterations  and  is  designed  to  force 
constraint  satisfaction  to  a predetermined  accuracy,  while 
the  norm  squared  of  the  variations  of  the  control  and  the 
parameter  is  minimized,  subject  to  the  linearized  constraints. 

The  principal  property  of  the  algorithms  is  that  they 
produce  a sequence  of  suboptimal  solutions,  each  satisfying 
the  constraints  to  the  same  predetermined  accuracy.  Therefore, 
the  values  of  the  functional  I corresponding  to  any  two  ele- 
ments of  the  sequence  are  comparable. 

The  gradient  phase  is  characterized  by  a descent  property 
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on  the  augmented  functional  J,  which  implies  a descent  pro- 
perty on  the  functional  I.  The  restoration  phase  is  charac- 
terized by  a descent  property  on  the  constraint  error  P.  The 
gradient  stepsize  and  the  restoration  stepsize  are  chosen 
such  that  the  restoration  phase  preserves  the  descent  property 
of  the  gradient  phase.  Hence,  the  value  of  the  functional  I 
at  the  end  of  any  complete  gradient-restoration  cycle  is 
smaller  than  the  value  of  the  same  functional  at  the  beginning 
of  that  cycle. 

Eight  numerical  examples  are  presented  to  illustrate 
the  performance  of  the  algorithms  associated  with  Problem  Pi 
and  Problem  P2.  The  numerical  results  show  the  feasibility 
as  well  as  the  convergence  characteristics  of  these  algorithms. 
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